Stress-driven phase transformation and the roughening of solid-solid interfaces 
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The application of stress to multiphase solid-liquid systems often results in morphological insta- 
bilities. Here we propose a solid-solid phase transformation model for roughening instability in the 
interface between two porous materials with different porosities under normal compression stresses. 
This instability is triggered by a finite jump in the free energy density across the interface, and it 
leads to the formation of finger-like structures aligned with the principal direction of compaction. 
The model is proposed as an explanation for the roughening of stylolites - irregular interfaces asso- 
ciated with the compaction of sedimentary rocks that fluctuate about a plane perpendicular to the 
principal direction of compaction. 

PACS numbers: 68.35.Ct, 68.35.Rh, 91.60.Hg 



Morphological instabilities in systems out of equilib- 
rium are central to most research on pattern formation. 
A host of processes give rise to such instabilities, and 
among the most intensively studied are the surface dif- 
fusion mediated Asaro-Tiller-Grinfeld instability Q}, Q, 
in the surfaces of stressed solids in contact with their 
melts, surface diffusion mediated thermal grooving and 
solidification controlled by thermal diffusion in the bulk 
melt [4]. In sedimentary rocks and other porous ma- 
terials local stress variations typically promote morpho- 
logical changes via dissolution in regions of high stress, 
transport through the fluid saturated pore space and pre- 
cipitation in regions of low stress. This phenomenon 
is known as pressure solution or chemical compaction. 
Such processes are often accompanied by the nucleation 
and growth of thin irregular sheets, interfaces or seams 
called stylolites jfj|. Stylolites form under a wide range 
of geological conditions as rough interfaces that fluctuate 
about a plane perpendicular to the axis of compression. 
They are common in a variety of rock types, including 
limestones, dolomites, sandstones and marbles, and they 
appear on scales ranging from the mineral grain scale 
to meters or greater. A common feature of stylolitic sur- 
faces is small scale roughness combined with large vertical 
steps in the direction of the compression. Residual unsol- 
uble minerals (i.e. clays, oxides) often accumulate at the 
interface as stylolites evolve. Despite the considerable at- 
tention given to the rich morphology of stylolites there is 
still no consensus on the mechanism(s) controlling their 
formation Q , , [i[ , @ • Here we demonstrate that even 
if the stylolite is a consequence of pressure solution alone, 
porosity or other material property gradients may drive 
the roughening process. In particular, we demonstrate 
that a compressional load normal to a solid-solid phase 
boundary gives rise to a morphological instability. 

Generally, rocks are heterogeneous bodies with spa- 
tially variable porosities. The strain energy densities may 
be larger in regions of high porosity (i.e. low modulus) 
than in regions of low porosity. Thermodynamically, the 
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FIG. 1: (Color online) Basic setup of the model for a mov- 
ing interface between two elastic solid phases, characterized 
by different Young's moduli (E2 > E\) and Poisson's ratios 
z/i , V2 ■ The interface boundary propagates with a normal ve- 
locity V n , when the solids are subjected to uniform far- field 
compressional stresses ctq in the vertical direction. 



total free energy of the system can be reduced by re- 
ducing the porosity variations. In this letter, a simple 
model for stylolite formation, in which high porosity rock 
is transformed into low porosity rock at the interface be- 
tween the low porosity and high porosity materials, is 
investigated. This solid-solid "phase transformation" is 
driven by the free energy density gradients, which can be 
substantial in regions with large porosity variations. The 
general approach used in this work could be applied to 
other solid-solid interface roughening phenomena. 

We consider a two-dimensional system divided into 
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elastic regions with different but homogeneous porosities 
(Fig. [I]). Without lack of generality, we limit our consid- 
eration to two dissimilar materials separated by a single 
interface. The stress boundary condition is a uniform 
compression in the vertical direction applied at the top 
and bottom boundaries. The two phases are separated 
by a sharp and coherent boundary, i.e. no defects or 
voids can form along the interface. This translates into 
continuity of the displacement vector u(r,t) across the 
interface, 



u 



0. 



(1) 



Here, and in other equations the brackets denote the 
jump in the quantity inside the brackets when the inter- 
face is approached from above and below. Under given 
load conditions the displacement field induced by the 
compression gives rise, in the linear regime, to a strain 
tensor of the form 
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The two solid phases are characterized by their Young's 
moduli, £a,2, and Poisson's ratios, v\^- When E\ < 
the upper region will be compressed the most and there- 
fore the elastic energy density will be higher in this re- 
gion. The first step towards a model for the roughening 
of a solid-solid interface is based on this simple obser- 
vation. First the elastic parameters of the materials are 
related to their porosities. Luo and Weng [10[ proposed 
a homogenization method relating the effective bulk and 
shear moduli to the porosity of the solid. Using this ap- 
proach, the effective Young's modulus decreases mono- 
tonically with the porosity. Consequently, a finite jump 
in porosity across the interface induces a jump in the 
elastic energy density, which drives the motion of the in- 
terface. Thermodynamically, the evolution of the inter- 
face corresponds to a phase transformation from a high 
to a low energy state. 

It is assumed that the phase transformation occurs on 
a time scale that is much longer that the time required 
for elastic waves to propagate across the system, and the 
system is therefore always in elastostatic equilibrium. For 
an isotropic and homogeneous elastic body the elastic 
equilibrium condition is given by 



dxj 



0, 



(3) 



together with the uniform uniaxial compression stress 
boundary condition. 



<Jij{x,y = oo ) = (JoSi y ySj y y < 0, 



(4) 



and the stress jump across the curved interface due to 
the effective surface tension is given by 



[cTijUj] = — jhzrii at the interface T^, 



(5) 



where k is the curvature and 7 is the local surface tension. 
In the limit of negligible surface tension, the stress vector 



is continuous across the interface ([cr^. 



0). 



For completeness, the basic principles used to derive 
an equation of motion are presented. When the system 
approaches an equilibrium configuration, the free energy 
will be a non-increasing function of time: 



d 
dt 
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Tds < 0, 
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where T is the free energy per unit volume and T is the 
interfacial free energy per unit area. Here, the subscript 
V indicates a volume integration and T t indicates inte- 
gration over the interface. The interfacial energy dissi- 
pation is obtained by confining the domain of integration 
to a narrow zone along the interface and taking the zero 
thickness limit [ll|. This gives 
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[T]V n ds j ( ^ - nTV n ) ds < 0, (7) 
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where V n is the normal velocity and n is the local cur- 
vature of the interface. This implies the differential form 
given by 



dt y 



[F])V n < 0, 



(8) 



where the first term is the total time derivative of the 
local interfacial energy density. The local free energy is a 
function of the surface tension only (like fluid-solid inter- 
faces) and thus, it is independent of time. Therefore, the 
time derivative can be neglected leading to the inequality 



(i&+[F\)V n <0. 



(9) 



In the linear response regime, this inequality is satisfied 
when the velocity (the thermodynamic flux) is a linear 
function of the driving force, (kJF+ [JF]), namely 



(10) 



with c > 0. In the absence of surface tension, the normal 
velocity is simply proportional to the jump in the strain 
energy, i.e. V n ~ [J 7 ]. While the dynamical law for the 
interface Eq. ([TO]) is very simple, the implementation in 
a numerical model is more challenging. 

The model (Fig. [T]) of the moving solid-solid inter- 
face was numerically implemented using the local force 
balance and energy dissipation equations (Eqs. © and 
([10])). The stress field is obtained using the Galerkin fi- 
nite element discretization of the elastostatic equations 
and the phase boundary is captured using the level set 
method. The level set method (|12j. |l3j ) is a power- 
ful and reliable technique for tracking surfaces in any 
number of dimensions. At any time t, the d-dimensional 
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FIG. 2: (Color online) Map of the logarithm of the elastic 
energy in the solids during the roughening process, with E\ — 
10 GPa, E 2 = 60 GPa, v x = v 2 = 0.3, and a = 0.05 MPa. 
Lower panel: initial h(x) at t = 0. Upper panel: interface at 
a later stage of roughening. 



interface T t may be defined as the zero level cut through 
a scalar field cp (d+l-dimensional surface), namely 



(/?(x, t) = , where x G Tt 



(ii) 



A change in the zero-level cut in response to a change 
in the scalar field may then be interpreted as a motion 
of the interface. Therefore the change in the scalar field 
must correspond to motion of the zero-level cut with a 
given normal velocity V, this is done by updating the 
scalar field using a simple advection-like equation 



~dt 



■V n \V<p\ = 0. 



(12) 



The advection of the level set function is solved on a sep- 
arate lattice using an upwind-scheme. The full dynami- 
cal model of the solid-solid phase transformation front is 
then given by this equation together with Eq. (fTU|) . 

The numerical simulations were started with a random 
interface generated by a directed random walk (Fig. [2 
lower panel). The temporal evolution of the interface 
(Fig. [21 upper panel) was then recorded for different ex- 
ternal stresses, <Jo, and elastic constants (E, v). The 
elastic constants were computed from homogenization re- 
lations between elasticity and porosity [l(j- 

Initially, the roughness of the interface grows exponen- 




FIG. 3: Roughness as a function of time for six different exter- 
nal compression stresses ao for a fixed jump in the Young's 
modulus and zero surface tension. The root means square 
height is plotted as a function of time, rescaled with the char- 
acteristic roughening time £*, on a semi-logarithmic scale. 
The data collapse shows the exponential roughening of the 
interface exp(£/£*) with a stress independent pre-exponential 
factor. Inset: the characteristic time as a function of external 
stress is given by t* ~ cr^ 2 . 
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h(s,t)-h(t)) ds ~ exp (£/£*), 



(13) 



with a characteristic roughening time £* that depends on 
the external stress and the jump in the elastic properties. 
To estimate the functional dependence of £* , a set of nu- 
merical simulations was performed. First, the external 
stress <Jo was systematically varied between 0.05 and 50 
MPa for fixed values of the elastic constants (E\ = 40 
GPa, E 2 = 60 GPa, v x = v 2 = 0.3). 

The results shown in Fig. [3] suggests that the charac- 
teristic time scales as £* ~ cr^~ 2 , and the prefactor de- 
pends on the elastic properties. In order to investigate 
this type of relation, the elastic constants across the in- 
terface were varied (E 2 = 50 GPa, E\ in the range 5-16 
GPa, v\ = v 2 = 0.3). The data for the interface roughen- 
ing collapses onto a curve which is exponential at small 
time-scales with a cross-over to a quadratic form at larger 
times (Fig. [4|). Extracting the characteristic time in the 
exponential growth regime and rescaling it with a 2 ,, the 
functional dependence with respect to the jump in the 
elastic constants across the interface was determined (see 
Fig. 31 inset). The cross-over from exponential to alge- 
braic roughening depends on the value of the jump. For 
large jumps in the elastic energy or porosity, the rough- 
ening quickly undergoes a transition from exponential to 
quadratic growth in time. This cross-over time is related 
to the formation and growth of the finger-like structures 
shown in Fig. [2l By a simple dimensional argument, we 
also determined that the generic dependence of £* on the 
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FIG. 4: Roughness as a function of time for different jumps in 
the Young's modulus or porosity [1/E] = 1/Ei — 1/E 2 ~ [4>] 
for six different external stresses and zero surface tension. The 
data is plotted against the time rescaled with t* estimated nu- 
merically from the initial exponential growth. At large times, 
the roughness grows quadratically with increasing time. In- 
set: Log- log plot of the characteristic time rescaled with the 
corresponding external stress ctq -t* as a function of the jump 
[1/E]. t* ~ erg • [1/E]~ a , with the scaling exponent a ~ 1.36. 



aligned with the direction of compression, and this may 
be explained by the model if the stress concentration at 
the fingers exceeds the yield strength of the material. 

To summarize, a simple solid-solid phase transforma- 
tion model that predicts a morphological instability of 
the interface under uniform compressional stress has been 
developed and investigated. The instability is triggered 
by a finite jump in the elastic properties across the inter- 
face and a concomitant jump in the free energy density. 
We also showed that the characteristic time of roughen- 
ing depends on the external applied stress and the elastic 
parameters jump, in such a way that a higher external 
compression load or a larger difference between the elas- 
tic properties of the phases shortens the time required to 
roughen the interface. This result allows the roughening 
time and formation rate of stylolites to be estimated as 
a function of burial depth in sedimentary basins. 

This project was funded by Physics of Geological Pro- 
cesses, a Center of Excellence at the University of Oslo. 
The authors are grateful to Paul Meakin for fruitful dis- 
cussions and comments. 



external stress and elastic constants, has the form 

*.„A„_L„* I JL (i4) 

where [J 7 ] = T\ — J~2 is the jump in the free energy density 
across the interface. The scaling relation is consistent 
with the numerical simulation results. The stress field 
may be calculated analytically with the same boundary 
and interface conditions for a straight interface perturbed 
by a sine function with small amplitude A. In linear per- 
turbative analysis (A « 1) without surface tension, the 
solution is obtained using the method of Airy's potentials 
for 2D elastostatics The energy jump [E] is propor- 
tional to kA 2 , where k = 2tt/X is the wave number. In 
other words, all the modes are unstable and those with 
the smaller wavelengths grow faster in the linear regime. 
The surface tension adds an ultraviolet cut-off resulting 
in small scales smoothening of the interface. The sys- 
tem was tested with and without surface tension and in 
both cases the qualitative behavior was the same - an ini- 
tial exponential roughening with a crossover to a finger- 
formation regime. Eventually, these fingers may stabilize 
due to transport of dissolved minerals and precipitation 
leading to pore clogging. Stress concentration at the tips 
is an important characteristic of the system. In stylolites 
the roughening is often accompanied by small fractures 
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